/**
 * @breif Construct matrix in hand
*/

label NC = mesh.nCells();

// 创建对应大小矩阵
simpleMatrix<scalar> A(NC);

// 为什么要自己初始化的原因：查看simpleMatrix构造函数的note
for(label i = 0; i < NC; i++)
{
    A.source()[i] = 0.0;
    for(label j = 0; j < NC; j++)
    {
        A[i][j] = 0.0;
    }
}

// 对角
for(label i = 0; i < NC; i++)
{
    A[i][i] = TEqn.diag()[i];
}


// 上下三角
for(label i = 0; i < TEqn.lduAddr().lowerAddr().size(); i++)
{
    // 1. 得到LDU地址
    label l = TEqn.lduAddr().lowerAddr()[i];
    label u = TEqn.lduAddr().upperAddr()[i];
    Info << "Address of lower: " << l << "\n upper: " << u << endl;

    // 2. 得到相应地址存储的值
    A[l][u] = TEqn.lower()[i];
    A[u][l] = TEqn.upper()[i];
}    

// 由边界带来的对矩阵的影响
    // 1. 遍历所有边界(下面两种都类似)
    // forAll(T.boundaryField(), patchI)
    forAll(mesh.boundary(), patchI)
    {
        // Info << "at BC---1: " << T.boundaryField()[patchI] << endl;   
        // Info << "at BC---2: " << mesh.boundary()[patchI].name() << endl;  
        const fvPatch &pp = T.boundaryField()[patchI].patch();
        forAll(pp, faceI)
        {
            label cellI = pp.faceCells()[faceI];
            A[cellI][cellI] += TEqn.internalCoeffs()[patchI][faceI];
            A.source()[cellI] += TEqn.boundaryCoeffs()[patchI][faceI];
        } 
    }


    Info << "\n==Coefficients of Matrix A==" << endl;
    for(label i=0; i<NC; i++)
    {
        for(label j=0; j<NC; j++)
        {
            Info<< A[i][j] << " ";
        }
            Info<< A.source()[i] << endl;
    }
    Info<< "\n==> Solution: " << A.solve() << endl;

        outputFilePtr.reset(new OFstream(outputDir/"matrix_dump_hand.txt"));

        outputFilePtr() << "Time: " << runTime.timeName() << endl;
        outputFilePtr() << "Matrix dump: " << endl;
            //打印及输出求解信息
        forAll(T,cellI)
        {
            outputFilePtr() << "(My construction)OF's solve() at cell: " << cellI << " : " << mesh.C()[cellI].component(vector::X)
                 << " with value: "<<T[cellI] << " "
                 << endl;
        }